###################################################################
###################################################################
############## Global Diffusion, Policy Flexibility, ############## 
################### and Inflation Targeting #######################
###################################################################
################## International Interactions #####################
###################################################################
######################### April 2019 ##############################
###################################################################
###################### Fabio Wasserfallen #########################
###################################################################
###################################################################


## Set working directory
setwd("/Users/.../ReplicationArchive/")


library(foreign)
library(arm)
library(texreg)
library(glmmML)
library(ggplot2)
library(pastecs)



################################################################
################# ANALYSIS OF THE ARTICLE ######################
################################################################


####################################################################
########### Figure 1: IT adoptions from 1990 to 2012 ###############
####################################################################

IT_adoption <- read.csv("IT_adoption.csv")

x <- IT_adoption$Year
x <- x[1:33]
y <- IT_adoption$Adoption
y <- y[1:33]
country <- as.vector(IT_adoption$Country)
country <- country[1:33]

par(mar=c(3,4.5,3,1)) 
plot(x,y,type="n", xlab="", ylab="Number of IT Adopters", main="Early Adopters                        Global Spread                                  ", cex.main= 1.15, font.main= 1.15, cex.axis=(1.1),cex.lab=(1.15), cex.sub=(1.1))
box(lwd=0.5)
text(x ,y ,labels= country, cex=(0.92))

abline(,,,1996.1,lty=2,lwd=0.5) 


############################################################
################# EMPIRICAL ANALYSIS  ######################
############ (Table 1, Figures 2 and 3) ####################
############################################################



########################################################################
##### Create ITflex Spatial Lag Variable with Inverse of Distance ######
########################################################################


### Read in Itflex data 
data.itflex <- read.csv("itflex_data.csv")
head(data.itflex)


### Prepare data for spatial lag creation
data.itflex$itflex[is.na(data.itflex$itflex)] <- 0
data.itflex$itrange[is.na(data.itflex$itrange)] <- 0
data.itflex$ittarget[is.na(data.itflex$ittarget)] <- 0
head(data.itflex)
length(data.itflex$itflex)



# Create yearly datasets and matrices of the ITflex, ITrange and ITtarget variable for all 27 years
data.1987 <- subset(data.itflex, year==1987)
data.1988 <- subset(data.itflex, year==1988)
data.1989 <- subset(data.itflex, year==1989)
data.1990 <- subset(data.itflex, year==1990)
data.1991 <- subset(data.itflex, year==1991)
data.1992 <- subset(data.itflex, year==1992)
data.1993 <- subset(data.itflex, year==1993)
data.1994 <- subset(data.itflex, year==1994)
data.1995 <- subset(data.itflex, year==1995)
data.1996 <- subset(data.itflex, year==1996)
data.1997 <- subset(data.itflex, year==1997)
data.1998 <- subset(data.itflex, year==1998)
data.1999 <- subset(data.itflex, year==1999)
data.2000 <- subset(data.itflex, year==2000)
data.2001 <- subset(data.itflex, year==2001)
data.2002 <- subset(data.itflex, year==2002)
data.2003 <- subset(data.itflex, year==2003)
data.2004 <- subset(data.itflex, year==2004)
data.2005 <- subset(data.itflex, year==2005)
data.2006 <- subset(data.itflex, year==2006)
data.2007 <- subset(data.itflex, year==2007)
data.2008 <- subset(data.itflex, year==2008)
data.2009 <- subset(data.itflex, year==2009)
data.2010 <- subset(data.itflex, year==2010)
data.2011 <- subset(data.itflex, year==2011)
data.2012 <- subset(data.itflex, year==2012)
data.2013 <- subset(data.itflex, year==2013)

itflex.matrix <- cbind(data.1987$itflex,data.1988$itflex,data.1989$itflex,data.1990$itflex,data.1991$itflex,data.1992$itflex,data.1993$itflex,data.1994$itflex,data.1995$itflex,data.1996$itflex,data.1997$itflex,data.1998$itflex,data.1999$itflex,data.2000$itflex,data.2001$itflex,data.2002$itflex,data.2003$itflex,data.2004$itflex,data.2005$itflex,data.2006$itflex,data.2007$itflex,data.2008$itflex,data.2009$itflex,data.2010$itflex,data.2011$itflex,data.2012$itflex,data.2013$itflex)

itrange.matrix <- cbind(data.1987$itrange,data.1988$itrange,data.1989$itrange,data.1990$itrange,data.1991$itrange,data.1992$itrange,data.1993$itrange,data.1994$itrange,data.1995$itrange,data.1996$itrange,data.1997$itrange,data.1998$itrange,data.1999$itrange,data.2000$itrange,data.2001$itrange,data.2002$itrange,data.2003$itrange,data.2004$itrange,data.2005$itrange,data.2006$itrange,data.2007$itrange,data.2008$itrange,data.2009$itrange,data.2010$itrange,data.2011$itrange,data.2012$itrange,data.2013$itrange)



### Create distance matrix (WDist)

library(cshapes)
wmat <- distmatrix(as.Date("2012-6-30"), type="capdist", useGW=F) 
#write.csv(wmat, file="wmat.csv") # done on February 6, 2019

# Take the 76 countries from the riskset.csv data (see Codes below)
vectorIDs <- c("339","160","371","900","305","211","145","571","140","355","20","402","155","100","94","344","316","390","130","92","366","375","220","255","452","350","90","91","310","395","750","850","205","666","325","51","740","501","367","368","212","343","820","70","565","790","210","920","93","475","385","770","150","135","840","290","235","360","365","345","317","349","560","732","230","780","380","225","800","52","640","200","2","165","101","551")

Wdist <- wmat[c(vectorIDs), c(vectorIDs)]
Wdist.inv <- 1/Wdist
diag(Wdist.inv) <- 0
nrow(Wdist.inv)
ncol(Wdist.inv)

# Row-standardize
W.dist.inv.RS <- matrix(NA,76,76)
for (i in 1:76) {
W.dist.inv.RS[i,] <- Wdist.inv[i,] / rowSums(Wdist.inv)[i]	
}


# Create the SLags

SLDISTitflex <- matrix(NA,76,27)

for (y in 1:27){
	for (i in 1:76){
		SLDISTitflex[i,y] <- W.dist.inv.RS[i,] %*% itflex.matrix[,y] 
		}
	} 
		
SLDISTitflex <- as.vector(t(SLDISTitflex)) 



SLDISTitrange <- matrix(NA,76,27)

for (y in 1:27){
	for (i in 1:76){
		SLDISTitrange[i,y] <- W.dist.inv.RS[i,] %*% itrange.matrix[,y] 
		}
	} 
		
SLDISTitrange <- as.vector(t(SLDISTitrange)) 




####################################################################
###############  Create dataset and run models #####################
####################################################################


### Attach SL variables to the riskset dataset for the analysis 
data.attach <- cbind(data.itflex, SLDISTitflex, SLDISTitrange)	

# Data for analysis
data <- read.csv("riskset_data.csv")
head(data)
data <- transform(data, t = data$year - 1986)
data <- transform(data, t_square = data$t^2/10, t_cubed = data$t^3/100, gdp_capita_in100k = gdp_capita/100000)
head(data)
summary(data)
	
	
# Merge 
data.analysis <- merge(data, data.attach, by=c("country_id","year")) # c() because many to one
head(data.analysis)

# Create index of IMF program dummies
data.analysis <- transform(data.analysis, imf_main = data.analysis$imf_sba + data.analysis$imf_eef + data.analysis$imf_saf + data.analysis$imf_prgf)



####################################################################
########### Table 1: Main Models of the Paper  #####################
####################################################################


mod1 <- glm(it_adoption.x ~ polity_2 + er_imf_fine + sumind_main_ext + pol_system + SLDISTitflex + fdi_in_pcgdp + exports_pcgdp + gdp_capita_in100k + imf_main + ka_open + t + t_square + t_cubed, family =binomial(link="logit"), data=data.analysis)
summary(mod1)

mod2 <- glm(it_adoption.x ~ polity_2 + er_imf_fine + sumind_main_ext + pol_system + SLDISTitflex + fdi_in_pcgdp + exports_pcgdp + gdp_capita_in100k + imf_main + ka_open + (polity_2*t) + t + t_square + t_cubed, family =binomial(link="logit"), data=data.analysis)
summary(mod2)

mod3 <- glm(it_adoption.x ~ polity_2 + er_imf_fine + sumind_main_ext + pol_system + SLDISTitflex + fdi_in_pcgdp + exports_pcgdp + gdp_capita_in100k + imf_main + ka_open + (polity_2* SLDISTitflex) + t + t_square + t_cubed, family =binomial(link="logit"), data=data.analysis)
summary(mod3)


screenreg(list(mod1,mod2,mod3), include.pvalues=TRUE ,stars = c(0.01,0.05,0.1), digits=3)
texreg(list(mod1,mod2,mod3), include.pvalues=TRUE ,stars = c(0.01,0.05,0.1), digits=3)



################################################################################
### Figure 2: Marginal Effect of Polity 2 score on IT Adoption over time ######
################################################################################

summary(mod2)
coef(mod2)
vcov(mod2)

# Marginal effect of polity if t = 1
Beta.POL.t1 <- coef(mod2)[2] + coef(mod2)[15]*1
SE.Beta.POL.t1 <- sqrt(vcov(mod2)[2,2] + 1^2*vcov(mod2)[15,15] + 1*2*vcov(mod2)[15,2])

# Marginal effect of polity if t = 13
Beta.POL.t13 <- coef(mod2)[2] + coef(mod2)[15]*13
SE.Beta.POL.t13 <- sqrt(vcov(mod2)[2,2] + 13^2*vcov(mod2)[15,15] + 13*2*vcov(mod2)[15,2])

# Marginal effect of polity if t = 27
Beta.POL.t13 <- coef(mod2)[2] + coef(mod2)[15]*27
SE.Beta.POL.t13 <- sqrt(vcov(mod2)[2,2] + 27^2*vcov(mod2)[15,15] + 27*2*vcov(mod2)[15,2])
     

t.score <- c(1:27)      
     
Beta.POL <- rep(NA,27)
SE.Beta.POL <- rep(NA,27) 
LOWER.CI.Beta.POL <- rep(NA,27)       
UPPER.CI.Beta.POL <- rep(NA,27)   

for(i in 1:27) {
Beta.POL[i]	<- coef(mod2)[2] + coef(mod2)[15]*i
}

for(i in 1:27) {
SE.Beta.POL[i]	<- sqrt(vcov(mod2)[2,2] + i^2*vcov(mod2)[15,15] + i*2*vcov(mod2)[15,2])
}

for(i in 1:27) {
LOWER.CI.Beta.POL[i] <- Beta.POL[i] - 1.96*SE.Beta.POL[i]
	}

for(i in 1:27) {
UPPER.CI.Beta.POL[i] <- Beta.POL[i] + 1.96*SE.Beta.POL[i]
	}


time <- c(1987:2013)

# Make Plot
layout(c(1,1),c(2,2), c(8,3.8))	
par(mar=c(6.2,4.3,0.5,2), mgp=c(3,0.8,0), cex.axis=1.2, cex.lab=1.2, font.main=1)
plot(time, Beta.POL, ylim=c(min(LOWER.CI.Beta.POL), max(UPPER.CI.Beta.POL)), xlab= "Year", ylab="Marginal Effect of Polity", col="black", cex.axis=1.1, cex.lab=1.15)
Axis(side=1, labels=F)
segments(time, LOWER.CI.Beta.POL, time, UPPER.CI.Beta.POL, col="black")
abline(0,0,lty=2)	



################################################################################################
### Figure 3: Marginal Effect of IT Flexibility SL on IT Adoption for Polity 2 score 0 to 10 ###
################################################################################################

summary(mod3)
coef(mod3)
vcov(mod3)

# Marginal effect of SL if Polity score = 0
Beta.SL.PS0 <- coef(mod3)[6] + coef(mod3)[15]*0
SE.Beta.SL.PS0 <- sqrt(vcov(mod3)[6,6] + 0^2*vcov(mod3)[15,15] + 0*2*vcov(mod3)[15,6])
Beta.SL.LOWER.CI <- Beta.SL.PS0 - 1.96*SE.Beta.SL.PS0
Beta.SL.UPPER.CI <- Beta.SL.PS0 + 1.96*SE.Beta.SL.PS0

# Marginal effect of SL if Polity score = 5
Beta.SL.PS5 <- coef(mod3)[6] + coef(mod3)[15]*5
SE.Beta.SL.PS5 <- sqrt(vcov(mod3)[6,6] + 5^2*vcov(mod3)[15,15] + 5*2*vcov(mod3)[15,6])
# Marginal effect of SL if Polity score = 10
Beta.SL.PS10 <- coef(mod3)[6] + coef(mod3)[15]*10
SE.Beta.SL.PS10 <- sqrt(vcov(mod3)[6,6] + 10^2*vcov(mod3)[15,15] + 10*2*vcov(mod3)[15,6])
     

Pol.Score <- c(0:10)      
     
Beta.SL <- rep(NA,11)
SE.Beta.SL <- rep(NA,11) 
LOWER.CI.Beta.SL <- rep(NA,11)       
UPPER.CI.Beta.SL <- rep(NA,11)   

for(i in 0:10) {
Beta.SL[i+1]	<- coef(mod3)[6] + coef(mod3)[15]*i
}

for(i in 0:10) {
SE.Beta.SL[i+1]	<- sqrt(vcov(mod3)[6,6] + i^2*vcov(mod3)[15,15] + i*2*vcov(mod3)[15,6])
}

for(i in 1:11) {
LOWER.CI.Beta.SL[i] <- Beta.SL[i] - 1.96*SE.Beta.SL[i]
	}

for(i in 1:11) {
UPPER.CI.Beta.SL[i] <- 	Beta.SL[i] + 1.96*SE.Beta.SL[i]
	}


summary(data.analysis$polity_2)	
d.polity_2 <- density(data.analysis$polity_2, na.rm=T, bw="nrd0", adjust=1)


# Make Plot
layout(c(1,2),c(2,2), c(8,3.8))	
par(mar=c(1.2,5.2,0.5,2), mgp=c(3,0.8,0), cex.axis=1.2, cex.lab=1.2, font.main=1)
plot(Pol.Score, Beta.SL, ylim=c(min(LOWER.CI.Beta.SL), max(UPPER.CI.Beta.SL)), xlab= "", ylab="", xaxt="n", col="black", cex.axis=1.1, cex.lab=1.15)
Axis(side=1, labels=F)
mtext(side=2, text="Marginal Effect of IT Flexibility", line=3.7, cex=1.15)
mtext(side=2, text="in Nearby Countries", line=2.5, cex=1.15)

segments(Pol.Score, LOWER.CI.Beta.SL, Pol.Score, UPPER.CI.Beta.SL, col="black")
abline(0,0,lty=2)	

par(mar=c(6,5.2,0,2), mgp=c(3,0.8,0), cex.axis=1.1, cex.lab=1.15, font.main=1)
plot(d.polity_2, xlim=c(0,10), xlab="Polity", ylab="Density", main="", cex.axis=1.1, cex.lab=1.15)
polygon(d.polity_2, col='grey')





##########################################################################
################# Analysis of the Online Appendix ########################
##########################################################################

##########################
##### Appendix A1 ########
##########################

# Summary Satistics
library(pastecs)
stat.desc(data.analysis$it_adoption.x)
stat.desc(data.analysis$SLDISTitflex)
stat.desc(data.analysis$SLDISTitrange)
stat.desc(data.analysis$polity_2)
stat.desc(data.analysis$er_imf_fine)
stat.desc(data.analysis$sumind_main_ext)
stat.desc(data.analysis$pol_system)
stat.desc(data.analysis$fdi_in_pcgdp)
stat.desc(data.analysis$exports_pcgdp)
stat.desc(data.analysis$gdp_capita_in100k)
stat.desc(data.analysis$imf_main)
stat.desc(data.analysis$ka_open)
stat.desc(data.analysis$cpi_movave_5y)
stat.desc(data.analysis$govparty)





##########################
##### Appendix A2 ########
##########################

# Create Box Plots for ITflex (Figure 1 in online appendix)

### Read in Itflex data 

data.itflex <- read.csv("itflex_data.csv")
head(data.itflex)


### Box Plots of IT flexibility for fully instututionalized and less democratic countries

ITFlexLessDemo <- data.itflex$itflex[data.itflex$polity2<6]
boxplot(ITFlexLessDemo)

ITFlexFullDemo <- data.itflex$itflex[data.itflex$polity2>9]
boxplot(ITFlexFullDemo)


# Organize dataset for this plot 
data.itflex.plot1 <- data.itflex[!is.na(data.itflex$polity2),]
data.itflex.plot1 <- data.itflex.plot1[!is.na(data.itflex.plot1$itflex),]
data.itflex.plot1 <- data.itflex.plot1[data.itflex.plot1$polity2<6 | data.itflex.plot1$polity2>9, ]
data.itflex.plot1$polity2
data.itflex.plot1$polity2[data.itflex.plot1$polity2<6] <- 0

# Make box plots (Figure 4 left plot)
data.itflex.plot1$polity2 <- as.factor(data.itflex.plot1$polity2)
p <- ggplot(data=data.itflex.plot1, aes(polity2, itflex)) +
geom_boxplot(outlier.size=1.5, outlier.shape=16, fill="gray") +
labs(title="", x="Polity Score", y="IT Flexibility")
p + scale_x_discrete(labels = c("0" = "< 6","10" ="10")) +
 theme_minimal() +
 theme(axis.text.x = element_text(size=12)) +
  theme(axis.text.y = element_text(size=12)) +
  scale_y_continuous(limits=c(2,16.5)) +
    theme(axis.title = element_text(size=14)) +
    theme(axis.title.x = element_text(margin = margin(18))) +
    theme(axis.title.y = element_text(margin = margin(0,18)))

median(data.itflex.plot1$itflex[data.itflex.plot1$polity2==0]) # 5.25
median(data.itflex.plot1$itflex[data.itflex.plot1$polity2==10]) # 4

mean(data.itflex.plot1$itflex[data.itflex.plot1$polity2==0]) # 6.03
mean(data.itflex.plot1$itflex[data.itflex.plot1$polity2==10]) # 3.97



# Organize dataset for this plot
data.itflex.plot2 <- data.itflex[!is.na(data.itflex$itflex),]
data.itflex.plot2$year[data.itflex.plot2$year>1989 & data.itflex.plot2$year<1997] <- 1 
data.itflex.plot2$year[data.itflex.plot2$year>1996 & data.itflex.plot2$year<2014] <- 2 

# Make box plots (Figure 4 right plot)
data.itflex.plot2$year <- as.factor(data.itflex.plot2$year)
p <- ggplot(data=data.itflex.plot2, aes(year, itflex)) +
geom_boxplot(outlier.size=1.5, outlier.shape=16, fill="gray") +
labs(title="", x="Stages of Diffusion", y="IT Flexibility")
p + scale_x_discrete(labels = c("1" = "1990-1996","2" ="1997-2013")) +
 theme_minimal() +
 theme(axis.text.x = element_text(size=12)) +
  theme(axis.text.y = element_text(size=12)) +
    scale_y_continuous(limits=c(2,16.5)) +
    theme(axis.title = element_text(size=14)) +
    theme(axis.title.x = element_text(margin = margin(18)))+
    theme(axis.title.y = element_text(margin = margin(0,18)))




####################################
##### Appendix A4, Table A4 ########
####################################

# Add Inflation and Governmnet Party as additional Controls to main models of paper (Table A4)

mod1.extended <- glm(it_adoption.x ~ polity_2 + er_imf_fine + sumind_main_ext + pol_system + SLDISTitflex + fdi_in_pcgdp + exports_pcgdp + gdp_capita_in100k + imf_main + ka_open + cpi_movave_5y + govparty + t + t_square + t_cubed, family =binomial(link="logit"), data=data.analysis)
summary(mod1.extended)

mod2.extended <- glm(it_adoption.x ~ polity_2 + er_imf_fine + sumind_main_ext + pol_system + SLDISTitflex + fdi_in_pcgdp + exports_pcgdp + gdp_capita_in100k + imf_main + ka_open + cpi_movave_5y + govparty  + (polity_2*t) + t + t_square + t_cubed, family =binomial(link="logit"), data=data.analysis)
summary(mod2.extended)

mod3.extended <- glm(it_adoption.x ~ polity_2 + er_imf_fine  + sumind_main_ext + pol_system + SLDISTitflex + fdi_in_pcgdp + exports_pcgdp + gdp_capita_in100k + imf_main + ka_open + cpi_movave_5y + govparty  + (polity_2* SLDISTitflex) + t + t_square + t_cubed, family =binomial(link="logit"), data=data.analysis)
summary(mod3.extended)


screenreg(list(mod1.extended,mod2.extended,mod3.extended), include.pvalues=TRUE ,stars = c(0.01,0.05,0.1), digits=3)
texreg(list(mod1.extended,mod2.extended,mod3.extended), include.pvalues=TRUE ,stars = c(0.01,0.05,0.1), digits=3)




####################################
##### Appendix A5, Table A5 ########
####################################

# Only Range as measure of IT Flexibility (Table A5)

mod1 <- glm(it_adoption.x ~ polity_2 + er_imf_fine + sumind_main_ext + pol_system + fdi_in_pcgdp + exports_pcgdp + gdp_capita_in100k + imf_main + ka_open + SLDISTitrange + (polity_2* SLDISTitrange) + t + t_square + t_cubed, family =binomial(link="logit"), data=data.analysis)
summary(mod1)

screenreg(list(mod1), include.pvalues=TRUE ,stars = c(0.01,0.05,0.1), digits=3)
texreg(list(mod1), include.pvalues=TRUE ,stars = c(0.01,0.05,0.1), digits=3)




##########################
##### Appendix A6 ########
##########################

# ER regime dummies (peg, intermediary, floating) instead of numerical variable 
# Create three ER regime dummies from er_imf_fine variable (three categories)
er_peg <- ifelse(data$er_imf_fine <6,1,0) # 1-5
er_intermediary <- ifelse(data$er_imf_fine >5 & data$er_imf_fine <11,1,0) # 6-10
er_floating <- ifelse(data$er_imf_fine >10,1,0) # 11-15

data.analysis <- cbind(data.analysis, er_peg, er_intermediary, er_floating)

# er_peg as reference category not included
mod3 <- glm(it_adoption.x ~ polity_2 + sumind_main_ext + pol_system + SLDISTitflex + fdi_in_pcgdp + exports_pcgdp + gdp_capita_in100k + imf_main + ka_open + er_intermediary + er_floating + (polity_2* SLDISTitflex) + t + t_square + t_cubed, family =binomial(link="logit"), data=data.analysis)
summary(mod3)


screenreg(list(mod3), include.pvalues=TRUE ,stars = c(0.01,0.05,0.1), digits=3)
texreg(list(mod3), include.pvalues=TRUE ,stars = c(0.01,0.05,0.1), digits=3)




